#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>
#include <vmath.h>
#include <xmltree.h>
#include <bvol.h>
#include <object.h>

#define VER_MAJOR	0
#define VER_MINOR	1

enum {
	ELEM_VERTEX,
	ELEM_NORMAL,
	ELEM_TANGENT,
	ELEM_TEXCOORD,
	ELEM_FACE,

	MAX_ELEM
};

static const char *elem_names[] = {
	"vertex",
	"normal",
	"tangent",
	"texcoord",
	"face",
	0
};

struct face {
	int vert_idx[4];
	int tc_idx[4];
	int norm_idx[4];
	int tang_idx[4];
};

struct mesh {
	vec3_t *vert;	/* vertex array */
	vec3_t *tc;		/* texture coordinate array */
	vec3_t *norm;	/* normal array */
	vec3_t *tang;	/* tangent array */

	int num_vert, num_tc, num_norm, num_tang;

	struct face *face;	/* face array */
	int num_face;

	struct bsphere bsph;
	struct baabox bbox;
};

struct tri_interp {
	vec2_t tex;
	vec3_t norm;
};

static void count_mesh_elem(struct xml_node *xml, struct mesh *m);
static void read_face(struct xml_node *xml, struct face *face);
static int elem_by_name(const char *name);

static inline int check_clock(vec3_t v0, vec3_t v1, vec3_t v2, vec3_t normal);
static int isect_triangle(ray_t ray, struct mesh *m, struct face *tri, struct isect_point *is);
static struct tri_interp interp_tri(struct mesh *m, struct face *tri, vec3_t ispos, vec3_t normal);

void plugin_info(FILE *fp)
{
	fprintf(fp, "mesh plugin (geometry) v%d.%d\n", VER_MAJOR, VER_MINOR);
	fprintf(fp, "  author: Ioannis Tsiombikas <nuclear@member.fsf.org>\n");
	fprintf(fp, "  description: adds support for tri/quad-mesh geometry.\n");
}

void free_data(void *data)
{
	struct mesh *m = data;

	if(m) {
		free(m->vert);
		free(m->tc);
		free(m->norm);
		free(m->tang);
		free(m->face);
		free(m);
	}
}

void *read_xml(struct xml_node *xml)
{
	struct mesh *mesh;
	struct xml_node *mesh_node, *c;
	int eidx[MAX_ELEM];

	mesh_node = xml->chld;
	while(mesh_node && strcmp(mesh_node->name, "mesh") != 0) {
		mesh_node = mesh_node->next;
	}
	if(!mesh_node) {
		return 0;
	}

	if(!(mesh = malloc(sizeof *mesh))) {
		goto err;
	}
	memset(mesh, 0, sizeof *mesh);

	count_mesh_elem(mesh_node, mesh);

	if(!(mesh->vert = malloc(mesh->num_vert * sizeof *mesh->vert))) {
		goto err;
	}
	if(!(mesh->tc = malloc(mesh->num_tc * sizeof *mesh->tc))) {
		goto err;
	}
	if(!(mesh->norm = malloc(mesh->num_norm * sizeof *mesh->norm))) {
		goto err;
	}
	if(!(mesh->tang = malloc(mesh->num_tang * sizeof *mesh->tang))) {
		goto err;
	}
	if(!(mesh->face = malloc(mesh->num_face * sizeof *mesh->face))) {
		goto err;
	}

	memset(eidx, 0, sizeof eidx);

	c = mesh_node->chld;
	while(c) {
		int elem;
		vec3_t val = {0, 0, 0};
		struct xml_attr *attr;

		elem = elem_by_name(c->name);

		if((attr = xml_get_attr(c, "id"))) {
			if(attr->type != ATYPE_INT) {
				fprintf(stderr, "mesh plugin: %s id field must be an integer\n", c->name);
				goto err;
			}
			eidx[elem] = attr->ival;
		}

		if((attr = xml_get_attr(c, "val"))) {
			if(attr->type != ATYPE_VEC) {
				fprintf(stderr, "mesh plugin: %s val field must be a vector\n", c->name);
			}
			val.x = attr->vval[0];
			val.y = attr->vval[1];
			val.z = attr->vval[2];
		}
		
		switch(elem) {
		case ELEM_VERTEX:
			mesh->vert[eidx[elem]++] = val;
			break;

		case ELEM_NORMAL:
			mesh->norm[eidx[elem]++] = val;
			break;

		case ELEM_TEXCOORD:
			mesh->tc[eidx[elem]++] = val;
			break;

		case ELEM_TANGENT:
			mesh->tang[eidx[elem]++] = val;
			break;

		case ELEM_FACE:
			read_face(c, &mesh->face[eidx[elem]++]);
			break;
		}
		c = c->next;
	}

	return mesh;

err:
	if(mesh) {
		free(mesh->vert);
		free(mesh->tc);
		free(mesh->norm);
		free(mesh->tang);
		free(mesh->face);
	}
	free(mesh);
	return 0;
}

struct xml_tree *write_xml(void *data)
{
	return 0;	/* TODO implement */
}

struct bsphere *get_bsph(void *data, mat4_t xform)
{
	return 0;	/* TODO implement */
}

struct baabox *get_bbox(void *data, mat4_t xform)
{
	return 0;	/* TODO implement */
}

void print_info(void *data, FILE *fp)
{
	struct mesh *m = data;

	fprintf(fp, "mesh: %d vertices, %d faces\n", m->num_vert, m->num_face);
}

int intersect(void *data, mat4_t xform, mat4_t inv_xform, ray_t ray, struct isect_point *isect_pt)
{
	struct mesh *m = data;
	int i;
	struct isect_point nearest, pt;

	/* bring the ray in object space */
	ray = ray_transform(ray, inv_xform);

	nearest.dist = DBL_MAX;

	for(i=0; i<m->num_face; i++) {
		if(isect_triangle(ray, m, m->face + i, &pt)) {
			if(pt.dist < nearest.dist) {
				nearest = pt;
			}
		}
	}

	if(nearest.dist != DBL_MAX) {
		if(isect_pt) {
			mat4_t inv_transp;
			m4_transpose(inv_transp, inv_xform);

			pt.pos = v3_transform(pt.pos, xform);
			pt.normal = v3_transform(pt.normal, inv_transp);

			*isect_pt = pt;
		}
		return 1;
	}
	return 0;
}

static int isect_triangle(ray_t ray, struct mesh *m, struct face *tri, struct isect_point *is)
{
	double n_dot_dir, t;
	vec3_t pos, normal, v1, v2, vo_vec;
	int i;

	v1 = v3_sub(m->vert[tri->vert_idx[1]], m->vert[tri->vert_idx[0]]);
	v2 = v3_sub(m->vert[tri->vert_idx[2]], m->vert[tri->vert_idx[0]]);
	normal = v3_cross(v1, v2);

	/* vertex-to-origin vector */
	vo_vec = v3_sub(ray.origin, m->vert[tri->vert_idx[0]]);
	
	if(fabs(n_dot_dir = v3_dot(normal, ray.dir)) < XSMALL_NUMBER) {
		return 0;	/* ray is parallel to the plane, no intersection */
	}

	/* calculate the intersection distance along the ray */
	t = -v3_dot(normal, vo_vec) / n_dot_dir;
	if(t < XSMALL_NUMBER) {
		return 0;	/* plane in the oposite subspace of the ray */
	}

	/* intersection point */
	pos.x = ray.origin.x + ray.dir.x * t;
	pos.y = ray.origin.y + ray.dir.y * t;
	pos.z = ray.origin.z + ray.dir.z * t;

	for(i=0; i<3; i++) {
		if(!check_clock(m->vert[tri->vert_idx[i]], m->vert[tri->vert_idx[(i + 1) % 3]], pos, normal)) {
			break;
		}
	}
	if(i < 3) {
		return 0;	/* the point is outside the triangle */
	}

	if(is) {
		struct tri_interp interp;

		is->dist = t;
		is->pos = pos;

		interp = interp_tri(m, tri, pos, normal);
		is->texcoords.x = interp.tex.x;
		is->texcoords.y = interp.tex.y;

		if(interp.norm.x == 0.0 && interp.norm.y == 0.0 && interp.norm.z == 0.0) {
			is->normal = normal;
		} else {
			is->normal = interp.norm;
		}

		if(v3_dot(is->normal, ray.dir) > 0.0) {
			is->normal.x = -is->normal.x;
			is->normal.y = -is->normal.y;
			is->normal.z = -is->normal.z;
		}
		is->normal = v3_normalize(is->normal);
	}
	return 1;
}

static inline int check_clock(vec3_t v0, vec3_t v1, vec3_t v2, vec3_t normal)
{
	vec3_t test = v3_cross(v3_sub(v1, v0), v3_sub(v2, v0));
	return v3_dot(test, normal) >= 0.0;
}


static void count_mesh_elem(struct xml_node *xml, struct mesh *m)
{
	struct xml_node *c;

	c = xml->chld;
	while(c) {
		int elem = elem_by_name(c->name);

		switch(elem) {
		case ELEM_VERTEX:
			m->num_vert++;
			break;
		case ELEM_NORMAL:
			m->num_norm++;
			break;
		case ELEM_TEXCOORD:
			m->num_tc++;
			break;
		case ELEM_TANGENT:
			m->num_tang++;
			break;
		case ELEM_FACE:
			m->num_face++;
			break;
		}

		c = c->next;
	}
}

static void read_face(struct xml_node *xml, struct face *face)
{
	int i = 0;
	struct xml_attr *attr;
	struct xml_node *face_chld = xml->chld;

	memset(face, 0, sizeof *face);

	while(face_chld) {
		if(strcmp(face_chld->name, "vref") == 0) {
			if((attr = xml_get_attr(face_chld, "vertex"))) {
				face->vert_idx[i] = attr->ival;
			}
			if((attr = xml_get_attr(face_chld, "normal"))) {
				face->norm_idx[i] = attr->ival;
			}
			if((attr = xml_get_attr(face_chld, "texcoord"))) {
				face->tc_idx[i] = attr->ival;
			}
			if((attr = xml_get_attr(face_chld, "tangent"))) {
				face->tang_idx[i] = attr->ival;
			}
			i++;
		}
		face_chld = face_chld->next;
	}
}

static int elem_by_name(const char *name)
{
	int i;
	for(i=0; elem_names[i]; i++) {
		if(strcmp(name, elem_names[i]) == 0) {
			return i;
		}
	}
	return -1;
}

#define LERP(a, b, t)	((a) + ((b) - (a)) * (t))

/* see Glassner; ch: 2; section: 3.3 (Eric Haines) */
/* TODO: readify this write-only piece of shit */
static struct tri_interp interp_tri(struct mesh *m, struct face *tri, vec3_t ispos, vec3_t normal)
{
	struct tri_interp interp;
	vec2_t res;
	vec3_t p00, p01, p11, p10, na, nb, nc, pa, pb, pc, pd;
	double du0, du1, du2, dv0, dv1, dv2;

	double s00, s10, s01, s11, s0, s1;
	double t00, t10, t01, t11, t0, t1;
	double nx00, nx10, nx01, nx11, nx0, nx1;
	double ny00, ny10, ny01, ny11, ny0, ny1;
	double nz00, nz10, nz01, nz11, nz0, nz1;
	/* TODO also do tangent */

	p00 = m->vert[tri->vert_idx[0]];
	p10 = m->vert[tri->vert_idx[1]];
	p01 = p11 = m->vert[tri->vert_idx[2]];

	s00 = m->tc[tri->tc_idx[0]].x; t00 = m->tc[tri->tc_idx[0]].y;
	s10 = m->tc[tri->tc_idx[1]].x; t10 = m->tc[tri->tc_idx[1]].y;
	s01 = s11 = m->tc[tri->tc_idx[2]].x; t01 = t11 = m->tc[tri->tc_idx[2]].y;

	nx00 = m->norm[tri->norm_idx[0]].x; ny00 = m->norm[tri->norm_idx[0]].y; nz00 = m->norm[tri->norm_idx[0]].z;
	nx10 = m->norm[tri->norm_idx[1]].x; ny10 = m->norm[tri->norm_idx[1]].y; nz10 = m->norm[tri->norm_idx[1]].z;
	nx01 = nx11 = m->norm[tri->norm_idx[2]].x; ny01 = ny11 = m->norm[tri->norm_idx[2]].y; nz01 = nz11 = m->norm[tri->norm_idx[2]].z;

	pa = v3_add(v3_sub(p00, p10), v3_sub(p11, p01));
	pb = v3_sub(p10, p00);
	pc = v3_sub(p01, p00);
	pd = p00;

	na = v3_cross(pa, normal);
	nb = v3_cross(pb, normal);
	nc = v3_cross(pc, normal);

	du0 = v3_dot(nc, pd);
	du1 = v3_dot(na, pd) + v3_dot(nc, pb);
	du2 = v3_dot(na, pb);

	dv0 = v3_dot(nb, pd);
	dv1 = v3_dot(na, pd) + v3_dot(nb, pc);
	dv2 = v3_dot(na, pc);

	if(fabs(du2) <= XSMALL_NUMBER) {
		double divisor = du1 - v3_dot(na, ispos);
		res.x = fabs(divisor) > XSMALL_NUMBER ? ((v3_dot(nc, ispos) - du0) / divisor) : 0.0;
	} else {
		vec3_t qux, quy;
		double dux, duy, ka, kb, u0, u1;
		
		qux = v3_scale(na, 1.0 / (2.0 * du2));
		quy = v3_scale(nc, -1.0 / du2);

		dux = -du1 / (2.0 * du2);
		duy = du0 / du2;

		ka = dux + v3_dot(qux, ispos);
		kb = duy + v3_dot(quy, ispos);

		u0 = ka - sqrt(ka * ka - kb);
		u1 = ka + sqrt(ka * ka - kb);

		res.x = (u0 < 0.0 || u0 > 1.0) ? u1 : u0;
	}

	if(fabs(dv2) <= XSMALL_NUMBER) {
		double divisor = dv1 - v3_dot(na, ispos);
		res.y = fabs(divisor) > XSMALL_NUMBER ? ((v3_dot(nb, ispos) - dv0) / divisor) : 0.0;
	} else {
		vec3_t qvx, qvy;
		double dvx, dvy, ka, kb, v0, v1;

		qvx = v3_scale(na, 1.0 / (2.0 * dv2));
		qvy = v3_scale(nb, -1.0 / dv2);

		dvx = -dv1 / (2.0 * dv2);
		dvy = dv0 / dv2;

		ka = dvx + v3_dot(qvx, ispos);
		kb = dvy + v3_dot(qvy, ispos);

		v0 = ka - sqrt(ka * ka - kb);
		v1 = ka + sqrt(ka * ka - kb);

		res.y = (v0 < 0.0 || v0 > 1.0) ? v1 : v0;
	}

	s0 = LERP(s00, s01, res.y);
	s1 = LERP(s10, s11, res.y);
	interp.tex.x = LERP(s0, s1, res.x);

	t0 = LERP(t00, t01, res.y);
	t1 = LERP(t10, t11, res.y);
	interp.tex.y = LERP(t0, t1, res.x);

	nx0 = LERP(nx00, nx01, res.y);
	nx1 = LERP(nx10, nx11, res.y);
	interp.norm.x = LERP(nx0, nx1, res.x);

	ny0 = LERP(ny00, ny01, res.y);
	ny1 = LERP(ny10, ny11, res.y);
	interp.norm.y = LERP(ny0, ny1, res.x);

	nz0 = LERP(nz00, nz01, res.y);
	nz1 = LERP(nz10, nz11, res.y);
	interp.norm.z = LERP(nz0, nz1, res.x);

	return interp;
}
